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ABSTRACT 

A design technique is proposed for linear regulators in which a feedback 
controller of fixed structure is chosen to minimize an integral quadratic 
objective function subject to the satisfaction of integral quadratic 
constraint functions. Application of a nonlinear programming algorithm to 
this mathematically tractable formulation results in an efficient and useful 
computer-aided design tool. Particular attention is paid to computational 
efficiency and various recommendations are made. Two design examples 
illustrate the flexibility of the approach and highlight the special insight 
afforded to the designer. 
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Introduction 


Research into control system optimization has followed a number of 
separate paths following the interest generated by the linear quadratic 
regulator (LQR) approach. Although it possesses attractive properties the LQR 
method suffers from two main drawbacks: 

i) the resulting controller requires access to the full set of plant 
states or to a state reconstructor, and, 
ii) the scalar quadratic cost function is inadequate for the 

representation of a set of design objectives. 

Levine and Athans (1970) pioneered the suboptimal regulator (SOR) 
approach by prespecifying the feedback controller structure but retaining the 
quadratic cost function adopted in the LQR design. Various features can be 
incorporated within the SOR context by expanding the cost function to include 
such measures as trajectory sensitivity and model-following (see Fleming 
1979). Computationally this is an attractive approach but it still suffers 
from the inadequacy of the scalar quadratic measure to describe all the 
different facets of system performance. 

Another line of attack was prompted by Schy, Adams and Johnson (1973) and 
Zakian (1973) in which system constraints and specifications are represented 
by a set of simultaneous algebraic inequalities. Thus the problem description 
could incude such basic control design parameters as overshoot, damping, 
settling time, etc. Polak and Mayne advanced this approach by posing semi- 
infinite programming problems in which design objectives are realized by 
minimizing a function subject to a set of inequalities where these objectives 
can be expressed as infinite dimensional constraints (Mayne, Polak and 
Sangiovanni— Vincentelli 1982) . Such problems require special mathematical 
programming algorithms to handle the infinite dimensional constraints and, 
sometimes, nondif ferentiable functions. 
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The linear quadratic constrained regulator (LQCR) method described here 
strikes a compromise between these two approaches. Common to all is the 
prespecification of controller structure; here an integral quadratic objective 
function is minimized subject to a set of integral quadratic constraint 
functions which may represent bounds on control energy, sensitivity measures, 
model-following errors, etc. An efficient solution procedure, based on 
readily available software, arises from this formulation and has led to an 
effective computer-aided design package being constructed on a 32-bit 
minicomputer. 

Two design examples which illustrate some novel features of the method 
are presented. Various controller configurations are studied for helicopter 
regulation where it is found that the identification of active constraints 
affords the designer additional insight into the problem. In a flight control 
example sensitivity reduction is an objective and a trade-off curve proves 
useful in selection of a suitable controller. 


2. A Typical LQCR Problem 

Given a linear time-variant plant 


where 


and 


x = A x + B u , 
-p p-p p-p 


y = c x 

— r» n — 


P-P 


X 

“P 


n x 1 plant vector 


u = m x l control vector 
— P 


y_^ - r x 1 output vector, 
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and a fixed controller structure, 


u 

-P 


K y 
o^p 




a typical LQCR design might seek a controller gain matrix, Kq, to regulate the 
output responses subject to certain control energy limitations, i.e. 

00 

T 

minimize / y Q n y dt, 
w.r.t.K 0 0" P °" P 

subject to 

oo 

2 

/ u dt < z , i=l , 2 , • • • ,m, 

0 * 

where z are constraint bounds. 

In a corresponding SOR design this problem would be approximately solved 
by successive minimizations of 


T* T 

J = / {yQy+uRu } d t , 

q l -p p^p -p p-p j 

where the designer strives to find the appropriate choice of Q p and R p to 
satisfy the control energy constraints. 

Thus what was previously solved in a number of unconstrained optimization 
designs (SOR) is now accomplished in a single constrained optimization design 
(LQCR). It is this "one-pass" solution procedure which makes the approach so 
appealing to a designer. Although an LQCR solution requires more computing 
time than an SOR solution it requires considerably less time than the sequence 
of SOR solutions necessary to solve the problem. Moreover use of the linear- 
quadratic formulation leads to suprisingly modest computing time overheads 
when compared with an SOR solution. 
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3. Design Procedure 

The design options available in the LQCR program are closely related to 
those found to be useful in an earlier SOR program (Fleming 1979). The design 
procedure summarized here is simply intended to be representative of the 
underlying concept. Future users will invent new design options appropriate 
to their needs. 

The linear time-invariant plant description is 


x = A x + B u % 
-P P-P p— P 


x (0) = x 
-P -P 


0 


y = c x 3 
-p v-p 


(1) 

( 2 ) 


and the controller configuration options available to the designer are: 
i) Full state feedback 


ii) Output feedback 


u = K x 
-P s-p 


u 

-P 

iii) Dynamic compensation 


K o^P 


U 

~P 


Ay + B x , 
cAp c— c ’ 


(3) 

(4) 

(5) 


where x 
equation 


is an s x l compensator vector which satisfies the dynamic 


x = C y + D x . 
— c c— p c— c 


x (0) = 0. 


( 6 ) 


Having selected a controller configuration the designer specifies the 
objective and constraint functions for the nonlinear programming problem: 
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such that 


minimize J r 


^ Z £ y ^ ^’^j***)*!* 


where z are constraint bounds and the set of cost functions, J , 

X X 

£=0,1, •••,q have the basic infinite-time quadratic integral form: 


J p = / (xJq_ x 41)^ u )dt, £=0,1, 

1 o ^ p r p ~ p p r p 


These cost functions may be expanded or modified according to the 

• • 

circumstances of the design. For example, a quadratic penalty term, x R x , 

c c £ -c 

may be included if the dynamic compensation option is selected. 

For the case where Ap and Bp (equation (1)) are functions of a scalar 
time-invariant parameter, a, a differential trajectory sensitivity 

vector, x g = 3x^/3a, is introduced which satisfies the equation 


x = Ax + Bu + A x + B 3u /3a, x (0) = 0, 
— s s— p s — p p— s p — p — s 


derived by partially differentiating equation (1) with respect to a, where 

A = 3A /3a and B = 3B /3a. If sensitivity reduction is an objective then 
s p s p 

T 

the quadratic sensitivity measure, x 0 x , may be included in the cost 

“® S £ S 

functions . 

In order to monitor how well a set of plant trajectories match a set of 
"model" trajectories either virtual model-following (VMF) or implicit model- 
following (IMF) terns may be Implemented. A quadratic measure of the 

T 

difference between plant and model outputs, (y ~y ) Q ( y -y ] , (VMF) or state 

P — ' ^ P ® 

, • • \ t r • 

derivatives, (x -x I Q lx -x I, (IMF) replaces the usual state measure, where 
p —nr p £ — p — m 

the model response is defined by 
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X 

— m 


A x 
m-m 


> 


y = c x , 

■=4n m-m 


x - v x 1 model state vector and y = r x 1 model output vector. 

Weighting matrices Q , Q , R and R will usually be diagonal and 

P SL H P £ i 

while it is possible to include more than one cost function option in either 
the objective or constraint functions it will be more beneficial in an LQCR 
design to identify these features separately. Indeed a cost function need 
only focus on a single vector element thus allowing the designer to exercise 
very precise control over competing functions in the minimization process. 


4. Formulation of Nonlinear Programming P roblem 

Depending on the selected design options the plant state vector, x will 

~P 

be augmented to incorporate a compensator state vector, x^, a trajectory 

sensitivity vector, x , and a model state vector, x to form the system state 

— s — m 

vector, x. This leads to the following concise expressions for the system 
state equation and cost functions: 

System x = A(K)x, x(0) = Xq, (7) 

where A(K) = A + + B 2 KC 2 (8) 

Objective/Constraint Functions 


-/ 

0 


x^OOxdt, 


£=0,1,* •• ,q , 


(9) 
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where 

Q £ (K) = q £ + C^K T R a KC,. (10) 

Matrices A,B^ ,B 2 ,C^ ,C 2 ,Q^ and R^ are easily derived from the input 

matrices A ,A ,A ,B ,B ,C ,C ,Q ,Q ,R and R ; their composition is 
P s’ p’ s’ p’ m P A s 4 * c £ P £ ’ 

described in Fleming (1979). Matrices A and are functions of the gain 

matrix, K, containing the optimization parameters, and it may have one of 
three constructions: 


i) Full state feedback 

ii) Output feedback 

iii) Dynamic compensation 

K = 

Note that it is not imperative for all of the elements, k^j, of K to be 
variable: some may be fixed to zero or constant values to aid investigation 

of gain redundancy, specific compensator structures, etc. (See Fleming 
1981). Collecting the variable elements of K into a parameter vector, k, we 
have the nonlinear programming problem: 



K = K 


K = K* 


minimize 
w.r.t. k 


( 11 ) 


subject to the inequality constraints 
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“ z £ < °> 1 = 1 , 2 ,»" ,q, 


( 12 ) 


where the cost functions are governed by the system state equation (7). 

Given a particular value for K the cost 

functions, J^,#. = 0,l,*»»q, (9), can be computed from 

J £ = tr ( p £ x 0 )» £=0,1, •••,q, (13) 

T 

where = XqXq anc * satisfies the Liapunov matrix equation 

P £ A + ATp £ = _Q £’ £=0,l,***,q, (14) 

provided that A is a stable matrix. 

Analytic expressions for the gradients of the cost functions (13) with 
respect to the controller gain matrix K, are derived here using an approach 
similar to that of Wilson (1970). Differentiating (13) with respect to a gain 
element, k^j, of matrix, K, we have 

3 J „ aT 3Q 

Tk^\ = 2tr ^ AP £,9k i ^ + tr ( A ’ (15) 

where A satisfies the Liapunov matrix equation 

AA T + AA + X Q = 0, (16) 


and combining (8), (10) and (15) it follows that 



2 ^ B 1 P £ AC 1 +B 2 P £ AC 2 +R £ KC 1 AC ? ‘ 


(17) 
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The gradient vector, 8J^/3k, can be easily constructed from (17). 

It is evident from (13) that the solution is initial condition dependent, 
however when initial conditions are unknown, matrix Xq may be modified so 

that either E{j }, ave J or max J. are evaluated (see Fleming 1979). 

y* x jL x 

- p 0 - p o 

It is a simple matter to extend the nonlinear programming problem, (11) 
and (12), to include direct constraints on gain parameters. Expressed as a 
general algebraic expression 

f (k) < 0, £=1 , 2 , • • • ,w (18) 

these constraints may simply represent bounds on certain variable gains, 
k^ j. Analytic gradients are obtained in a straightforward manner. 


5. Solution Procedure 
5.1. Nonlinear Programming Algorithm 

Following a short survey of algorithms using the ADS program package 
(Vanderplaats 1983) it was established that the use of analytic gradients (17) 
is preferred to that of using finite difference approximations, resulting in 
improved accuracy and solution times. Four algorithms were compared: 

a) Method of feasible directions (Zoutendijk 1960; Vanderplaats 1973), 

b) Sequential unconstrained minimization technique (SUMT) using the 
quadratic exterior penalty function method (Fiacco and McCormick 1968) , 

c) SUMT using the quadratic extended interior penalty function method 
(Haftka and Starnes 1976), and 

d) Augmented Lagrangian multiplier method (Powell 1978). 
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Although (c) was the most accurate, (a) was found to be the most satisfactory 
algorithm within the CAD context, yielding solution times 6-8 times faster 
than (b), (c) and (d) for comparable accuracy levels. The promising iterative 
quadratic programming method of Powell (1978) will be tested using a 
forthcoming version of ADS. 

Zoutendijk's Method of Feasible Directions initially finds a feasible 
point and then proceeds by iteratively searching along feasible directions. 
In the program package the BFGS variable metric method (e.g. Fletcher 1970) is 
employed within the feasible directions method if no constraints are 
violated. At each iteration gradient information is required only for the 
active or violated constraints. The objective and constraint functions are 
computed at each iteration and within the line search procedures for each new 
estimate of K. 


5.2 . Liapunov Matrix Equations (LMEs ) 

A breakdown of CPU usage reveals that function and gradient evaluations 
dominate the solution time and these evaluations, in turn, are subject to the 
efficiency of the solution of the LMEs, (14) and (16). Generally recognized 
as the most efficient general LME solver, the method of Bartels and Stewart 
(1972) possesses additional properties which may be exploited in this 
application. It is a transformation method which reduces A, (14), to its 

real Schur form and then obtains a solution by solving sets of linear systems 
whose individual orders do not exceed four. Therefore only one Schur 
reduction of A per estimate of K is required for the solution of the 
(q+2) LMEs, (14) and (16), when computing the objective and constraint 
functions, (13), and their gradients, (17). The solution of the LMEs is 
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further accelerated by exploiting certain structural and sparseness properties 
of A,Q^ and Xq. The computing time for the Bartels-Stewart algorithm is 

Q 

proportional to N , where N is the order of the LME. 

It should be noted that there are circumstances under which the LMEs are 
more efficiently solved using the direct method of MacFarlane (1963). 
Although the computing time for this method for one LME solution is 
proportional to N^, subsequent solutions for different (14) rely only on 

the application of back substitution. However the execution time superiority 
of this method only prevails for certain combinations of small N and large 
q (12). For large N the Bartels-Stewart algorithm is the most efficient 

and is the recommended technique. 


5.3. Obtaining and Maintaining a Stable Matrix A 

While it is permissible to start the optimization from an infeasible 
point with respect to the constraints (12) and (18), the solution of (14) 
demands that the initial parameter vector, k, and subsequent estimates of this 
vector stabilize the system (7). An expression constraining eigenvalues of 
A to have negative real parts cannot be included in the nonlinear programming 
problem description since violation of this constraint invalidates the 
solution of (14). A steepest descent technique similar to that of Koenigsberg 
and Frederick (1970) is recommended to search for a stabilizing value of k 
should one be unavailable. 

Subsequently the eigenvalues of A are monitored throughout the 

optimization search procedure. Although the problem is such that the routine 
will tend to generate stabilizing values of _k, computational traps have been 
set to inhibit excursions into the unstable region. Should such a violation 
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occur the line search step is repeatedly halved until a stable value of k is 
reached. Use of the Bartels-Stewart method for solving LMEs permits economic 
evaluation of the eigenvalues from the diagonal and principal subdiagonal 
elements of the Schur form of A. 


6. Application Examples 

6.1 Alternative Controller Configurations 

Helicopter longitudinal dynamics are described in Michael and Farrar 

(1973) for a plant which has four states, Xp = [y x ,y z ,0 , 0 ] T , and two 

controls, u = [u 1 ,u ? ] T . The goal is to satisfactorily regulate y and y 
P 1 z x z 

(forward and vertical velocities) employing feedback from only two states, 

0 and 0 (pitch angle and pitch rate) , since the use of airspeed sensors for 
anc * un desirable. An LQR design, incorporating full state 

feedback, defines a satisfactory model response: 


where 


x = (a +B K*) 

— m v p p c' 


— m 


X 

s y — m 


= A x 


* 

u = K x . 
-p s-p 


The aim here is to attempt to match this model response without recourse to 
the use of airspeed sensors and recognizing control magnitude limitations. 

Thus we have the LQCR design problem in which we seek a controller having 
an acceptable configuration to minimize the model-following error term 

00 

J n = / (x -x ) T Q (x -x ldt, 

o o — p m; pV— p -In ' 
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where Q p = diag (1,1, 0,0), subject to the control energy constraints 


J - = / u?dt - 4.06 < 0, 

1 0 

00 2 

J = / u.dt - 6.57 < 0, 

0 

where the control energy bounds are those limiting the LQR design. Tabulated 
design results for three controller configuration besides the LQR controller 
are presented. 

Table 1 . Design results for various controller configurations 


Controller 


LQR u 

-P 


1 u 

-P 


Gain Matrix, K 


-2.56 

0.38 


0.29 0.50 0.24 

1.98 -0.26 -0.06 


x 

-P 


0 


0.16 

-0.25 


0.14 

0.22 



0.072 


Active Constraints 


J l> J 2 


J 2 


2 


u 

-P 


0.18 

-0.25 


9 


0.35 


J 


2 


u 


" 0.32 -0.24 


"e 

-p 

= 

-0.06 -0.30 


• 

X 

c I 


1 -1.37 


X 

c 


Identification of the active constraint for the design of Controllers 1 and 2 
indicates that relaxation of the control energy bound on U 2 will lead to 
improved model— following. It is clear from Figs. 1 and 2 that excellent 
model-following is achieved by Controller 1, while Controller 2, employing 
pitch angle feedback only, obtains less satisfactory results. 
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Implementing controller dynamics with pitch angle feedback (Controller 3) 
produces model-following results comparable to those of Controller 1. For 
this design gain element k 31 was fixed at unity in order to realize a 
minimal parameter form for the compensator. 

It is interesting to find that for the control configurations 1 and 2, 
only control u 2 is exercised to the limit (its energy term is active at the 
minimum) in pursuit of the model-following objective, thus implying that 
increased u-^ control energy degrades the model-following capability. This 
is simply caused by implementation of these new control configurations and is 
not due to any properties of the model since the same effect was noted for an 
objective function containing plant output terms alone. 

A fuller discussion of this example is contained in Fleming (1983), where 
it is shown that the control energy bounds are effective in limiting maximum 
control magnitudes. 


6.2 Sensitivity Reduction 

This example exercises two cost function options: sensitivity reduction 

and implicit model-following. IMF is employed in favor of VMF since it 
results in a lower-order system description and is sufficiently accurate for 
this application. 

We consider the flight dynamics of an aerodynamically unstable aircraft 
(Grubel and Kreisselmeier 1974): 

^ = y a )*p + y«)u, 


where the components of x^ are incremental longitudinal velocity, flight 

path angle, y, pitch rate and pitch angle. The control input is the elevator 

deflection. Plant matrices Ap and B p are dependent on a parameter 

a (0.3<a<0.7) which is the relative position of the c.g. of the aircraft. 

Plant and sensitivity matrices Ap, B p , A g , and B g are evaluated for 

dp = 0.5 and may be determined from Grubel and Kreisselmeier (1974). In 

their paper, the authors solved the LQR problem and observed the flight path 

T 

angle response to the initial condition x =[0101] (see Fig. 3). The 

p o 

nominal response (a^ = 0.5) satisfies the design goal of following a step 
input with essentially no overshoot in a 5 percent settling time of 2s. 

However the trajectories for off-nominal parameter values are 

unsatisfactory: a sluggish response for a = 0.3 and 73 percent overshoot 

for a = 0.7. 

The aim of the LQCR design is to find a controller of similar magnitude 

which gives a comparable flight path angle response for a = 0.5 and is less 

sensitive to parameter variations. We therefore seek the full state feedback 

controller, u = Kx , which minimizes the sensitivity measure 
-P ~P 


T 

J n = / xQxdt 
0 J 0 — s x s— s 


(19) 


subject to the constraints on model-following errors, 


J. = / fx -x )q (x -x Idt - z. < 0 
1 J p — m ; p v — p -nr 1 


( 20 ) 


and control energy, 
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where Q g = diag (0, 1,0,0), Q p = I,R p = 100, z 2 = 7.77 and the model state 
vector is derived from the LQR response for a = 0.5. The bound z 2 
corresponds to the LQR control energy measure and the model-following error 
bound, z 1# is open to experimentation by the designer. It was found that the 
control energy constraint was ineffective in limiting the magnitude of 
u p (0) which played an important role in sensitivity reduction. To maintain a 
consistent comparison with the LQR design the following two algebraic 
constraints of the form (18) were added: 


k 12 + 


14 


- u 


< 0 


( 22 ) 


where 


( k 12+ k i 4 ) 


u < 

p o 


0, 


(23) 


u (0) = K x 
P -P 


0 



0 

1 

0 

1 


= k 12 + k 


14’ 


and for the LQR design, u (0) = u = 0.60. 

P P 0 

LQCR designs were carried out for a range of values of z 1 to examine 
the expected trade-off between sensitivity reduction and model-following (see 
Fig. 4). This trade-off curve suggests that Controller A is a candidate for 
"best" design and the corresponding flight path angle response is illustrated 
in Fig. 5(a). Here trajectory dispersion of y is reduced with overshoot for 
a — 0.7 less than 35% while the nominal response has no overshoot and a 5 
percent settling time of 3.0s. Controller B improves on model-following 
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capability at the expense of sensitivity reduction (see Fig. 5(b)). In both 
designs the active constraints are equations (20) and (22). 

At this level of model-following error (z^l) we investigate the effect 
of abandoning the control constraints altogether, i.e. dropping constraints 
(21)-(23) , to find the amount of sensitivity reduction possible under these 
relaxed conditions. The resulting flight path angle y, response due to 
Controller C is illustrated in Fig. 5(c), where we observe, in particular, 
that overshoot for the a = 0.7 case is reduced to 44%; the actual 
sensitivity measure (19), Jq, is 3.40. Since the additional control effort, 

00 

/ 100u 2 dt = 8.09, 

0 P 

u (0) = 0.76, 

P 

is relatively minor the implication is that relaxing control constraints pays 
substantial dividends for this sensitivity-accented problem. 


6.3. Computing Details 

These examples were worked on a VAX 11/95 minicomputer which operated a 
conversational-mode CAD package employing the algorithms recommended in 
Section 5. Both examples required augmented system descriptions (7) of order 
8 (order 9 for the dynamic compensation case in 6.1) and the number of 
minimizing variables ranged from 2 to 5. Computing times for individual 
designs varied according to solution accuracy demands and initial controller 
gain estimates but took typically 20-100 seconds of CPU time. Tests carried 
out on corresponding S0R designs indicated that LQCR designs required 4-8 
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times more computing time. It must be remembered, however, that the SOR 
design process is essentially iterative requiring a large number of designs 
before approaching the results of a single LQCR design. 


7 . Concluding Remarks 

Nonlinear programming has been applied to regulator design in a variant 
of the LQR design procedure. The approach has significant advantages allowing 
different controller configurations to be tested and invoking sensitivity and 
model-following terms together with the more usual state and control terms in 
the design functions. It achieves design goals in a more direct and 
convenient manner than its SOR design counterparts. Used as a one-pass 
solution procedure (Section 6.1) or in the generation of trade-off curves 
(Section 6.2) each design yields a variety of information from inspection of 
objective and constraint values at the solution. 

Its format is flexible, based on a linear— quadratic formulation, and may 
be easily modified for different user's requirements, e.g. inclusion of 
disturbance measures, new sensitivity approaches. The integral quadratic 
measures which represent the objective and constraint functions may be readily 
interpreted as RMS values in designs which include a noise term in the system 
description. 
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Flight path angle, y , responses for 
b) Controller B, c) Controller C. 


a) Controller A 
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